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SUBMERGED  STRUCTURAL  RESPONSE  TO  WEAK  SHOCK  BY 
COUPLED  THREE-DIMENSIONAL  RETARTED  POTENTIAL  FLUID 
ANALYSIS-FINITE  ELEMENT  STRUCTURAL  ANALYSIS 


INTRODUCTION 

This  study  is  s  continuation  of  the  effort  to  apply  the  retarded  potential  integral  to  increasingly 
general  classes  of  problems  in  Quid-structure  interaction.  Mitzner  (1.2]  developed  a  method  for 
discretizing  the  integral  spatially  and  temporally,  and  applied  the  technique  to  the  problem  of  a  spheri¬ 
cal  rigid,  immovable  body  impinged  by  a  plane  pressure  wave  with  a  smooth  waveform.  Huang, 
Everstine,  and  Wang  P]  extended  the  method  to  handle  a  non-rigid  boundary  condition,  and  tested 
the  case  of  a  spherical  elastic  shell,  modeled  with  axisymmetric  shell  elements,  impinged  by  a  plane 
saep  wave.  However,  the  special  problem  posed  for  the  integral  formulation  by  a  pressure  discon¬ 
tinuity  was  by  solving  the  rigid  scattering  portion  of  the  problem  closed-form  using  separation 

of  variables.  Neilson,  Lu,  and  Wang  (7]  directly  tackled  the  problem  of  the  discontinuous  shock 
wave  front  for  rigid  body  problems  and  axisymmetric  surfaces  using  the  wavefront  line  integral  term. 

This  work  builds  upon  the  computer  program  developed  by  Mitzner  and  extended  to  non-rigid 
boundaries  by  Huang,  Everstins,  and  Wang.  A  fully  thrse-dimsnsional  boundary  element  formulation 
for  arbitrary  shaped  surfaces  has  bean  implemented.  Tbs  entire  method  is  designed  to  be  com¬ 
patible  with  finite  element  methods  and  has  been  linked  to  the  ADINA  Finite  element  program.  Con¬ 
tinuous  pressure  loadings  can  be  successfully  used  u>  approximate  loading  discontinuities,  avoiding  the 
difficulties  inherent  in  the  wavefront  line  integral  approach. 

A  concurrent  effort  in  the  creation  of  a  general  three-dimensional  analysis  capability  is  the 
development  of  a  parametric  patch  surface  geometry  definition  for  the  retarded  potential  integral.  This 
allows  arbitrarily  shaped  surfaces  to  be  specified  by  the  coordinates  of  a  finite  number  of  points,  as  is 
done  with  finite  elements.  This  development  is  described  in  a  related  report  [10],  and  is  used  in  the 
computations  described  in  this  report. 

RETARDED  POTENTIAL  FORMULATION 


The  retarded  potential  integral  is  the  solution  of  the  linear  wave  equation. 


1  &p 

c2  dt2 


subject  to  the  boundary  condition,  dp  /dn  *■  -  p  w .  This  is  the  governing  differential  equation  for 
a  linearly  compressible,  inviscid  fluid  with  negligible  changes  in  density.  The  assumption  of  linear 
fluid  behavior  is  said  to  be  accurate  for  pressures  up  to  30,000  pst  [9]. 

A  discussion  of  the  retarded  potential  equation  and  its  discretized  approximation  as  applied  to  the 
submerged  structural  problem  may  be  found  in  Refs.  (1,2.3).  A  good  derivation  is  found  in  Ref.  (4). 
The  form  of  the  equation  for  calculation  of  total  pressures,  p,  on  a  bounding  surface  in  an  infi¬ 
nite  fluid  subject  to  a  continuous  incident  pressure  Field  loading,  .  is  as  follows: 
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ppj)  -  2 p*&J)  ~  ptlr  fs  dS' 

{'<*••''> +  tI^™}  f£"' 

where  t'  •  t  -  — 
c 

and  R  -  \I  -  r| 


5  is  the  position  of  a  field  point  on  the  surface, 

J'  is  the  location  of  an  integration  point  on  the  surface, 

r'  is  the  retarded  time  of  the  influence  of  field  quantities  and  the  integration  point, 

p  is  the  fluid  density, 

e  is  the  sonic  velocity  in  the  fluid, 

m  '  represent!  the  surface  normal  direction  at  an  integration  point  oriented  out  of  the  fluid, 
w'  is  the  normal  direction  acceleration  of  the  surface  oriented  out  of  the  fluid. 

IMPLEMENTATION 

For  computation,  the  retarded  potential  integral  must  be  discretized  and.  for  a  non-rigid  bound¬ 
ing  structure,  linked  to  a  structural  analysis  code.  The  surface  pressure  field  is  approximated  by  sub¬ 
dividing  the  surface  into  zones  of  constant  pressure  coinciding  with  individual  finite  element  surfaces, 
and  with  the  surface  normal  acceleration  field  also  approximated  as  constant  on  the  same  zonal  sur- 
thce.  This  normal  acceleration  is  obtained  from  the  finite  element  computation  by  averaging  (he 
acceleration  vectors  from  the  comer  nodes  of  each  finite  element  and  using  the  component  normal  to 
the  censer  of  the  zone.  Pressures  are  discretized  in  time  to  coincide  with  the  time  steps  of  the  finite 
element  program.  Prior  to  conducting  the  time  history  computation  of  pressures  and  structural  vibra¬ 
tions,  a  matrix  of  influence  coefficients  must  be  calculated  that  expresses  the  dependence  relations  of 
current  pressures  on  past  pressures  and  past  and  present  accelerations. 

First,  the  time  derivative  of  pressure  must  be  formulated  in  terms  of  pressures.  A  three-point 
backward  difference  formula  is  used  in  this  implementation.  Then,  for  each  zone,  pressure  and 
acceleration  can  be  factored  out  of  the  integral  expressions.  There  is  a  need,  however,  to  subdivide 
the  zones  into  a  mesh  of  subzone  elements  in  order  to  accurately  numerically  integrate  the  geometric 
influence  factors,  as  well  as  to  accurately  obtain  the  time-delay  of  pressure  influences  from  all  pans 
of  the  structure.  For  the  purpose  of  measuring  the  distance  between  integration  and  field  points  and 
thereby  the  tune-retardation,  the  field  point  is  located  at  the  center  of  the  field  zone  and  the  integra¬ 
tion  point  is  located  at  the  center  of  the  subzone  integration  element.  The  exception  to  this  rule  is 
that  the  influence  of  the  singular  integration  element,  which  surrounds  the  field  point,  is  lumped  ai  a 
distance  from  the  field  point  of  one-fourth  the  diagonal  dimension  of  the  singular  element.  The 
integration  points  will  not  fail  at  a  whole  number  of  time  intervals  from  the  field  point,  so  an  interpo¬ 
lation  from  adjacent  whole  time  intervals  is  needed  to  determine  the  pressure  at  the  integration  point 


Fundamental  Relationship*  in  the  Discretization 


unman*. 

X  -  petition  vector 

R  -  distance  between  integration  and  field  points 
j  -  field  point  (and  zone)  index 
k  -  insagration  zone  index 
l  -  subcone  element  index 
(kJ)  -  integration  point  indices 
c  -  wav*  velocity 
t  •  enrrant  time 
f '  -  retarded  time 
r  •  time  ssep  size 
m  -  time  step  index 

m  ’  •  time  step  index  of  retarded  time  (not  a  whole  number) 
i  -  generic  index  of  time  retardation  or  delay 
/  'jh  -  number  of  time  steps  of  delay  in  the  influence  of  integration  point  (k,l) 
on  field  point  j  (not  a  whole  number) 


Pressure: 

Pjm  -P<Xjjnr) 

Ptm  •  P<XmJ*'t) 


Acceleration: 

"jm  •  »(f/i*r) 
wfe,-  ■  w(T*.nrr) 


Rjh  -  |  Xj  -Iu  | 

Retarded  Tune: 

t'  •  t  -  R/c 

m'r  ■air  -  fl^/c  -  air  -  i '  ^  r 

Tims  Retarded  Pressure,  Acceleration,  and  Pressure  Derivative  in  terms  of  pressures  and  acceleration 

at  whole  time  steps: 
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wham,  for  linaar  interpolation  of  field  quantities  between  whole  time  steps, 
fpM  “  1  ~  1  >  i-i'ju  >  -1 

fjut  ■  0,  otherwise 

sad,  for  the  3-point  backward  difference  discretization  of  the  pressure  derivative, 
f}i  ■  \fm  ~  Vjnu- 1  +  y  fjuj-i 


Retarded  Potential  Integral: 
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Inserting  the  above  expressions  for  time  retarded  pressures,  accelerations,  and  pressure  derivatives 
into  the  discrete  Retarded  Potential  Integral,  summing  over  the  l  index,  and  consolidating  the  time- 
independent  factors  into  matrices  called  A  and  B,  results  in  an  expression  of  the  form, 

0  Pjm  m  2Pjm  +  2-  "kjm-i  +  £  &jU  Pkjm-i 

kj  kj 


which  allows  the  current  pressures  pJm 


to  be  computed  from  past  pressures  and  current  and  past 


The  particular  approximations  used  for  the  geometric  influence  factors. 


ft  -LJ*dr 

LK  &  tot'  " 


k  R 


and 


ju 


,  depend  on  whether  or  tax  the  integration  element,  (£./),  is  singular  with 


to  field  point,  j .  This  is  discussed  further  in  Ref.  10. 


The  numerical  integration  over  the  subzone  elements  is  a  subject  left  to  the  related  paper  con¬ 
cerning  a  parametric  surface  patch  basis  for  the  geometric  calculations  of  the  method  [10].  The  origi¬ 
nal  formulation  of  the  numerical  integration  is  presented  in  [1,2]. 


A  staggered  solution  technique  is  used  in  the  analysis  so  that  for  each  time  step  the  structural 
and  fluid  dynamic  problems  can  be  solved  alternately  rather  than  concurrently,  greatly  simplifying 
computation,  yet  nearly  preserving  the  simultaneity  of  the  two  problem  solutions.  For  the  step  by 
step  soiunou  of  the  fluid-structure  interaction  problem,  the  cycle  of  computation  is  as  follows: 


0.  At  tune  step  0,  initialize  pJ0  •  0 


1.  At  time  step  m,  tm  m  mr.  Apply  pressures,  pjm  to  structural  model. 

2.  Calculate  accelerations,  at  time  from  structural  model. 

3.  Calculate  pressures,  pj  m  +  t,  at  time,  tm  +  {,  with  integral  equation,  due  to 

(the  time  retarded  surface  accelerations;  and  (the  time  retarded  surface  pres¬ 

sures). 


4.  If  m  + 1  is  the  last  desired  time  step,  stop  the  solution. 
3.  Increment  m.  Go  to  computational  step  l. 
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Tbs  computational  cycle  is  depicted  in  Fig.  1,  showing  the  relationships  among  the  processes 
sod  the  date  input  and  generated.  From  the  point  of  view  of  the  retarded  potential  calculations  there 
an  two  iwHapmHui*  inputs:  the  incident  wave  loading  and  the  surface  accelerations.  It  is  sometimes 
co  think  separately  of  the  pressure  field  contributions  of  these  inputs.  For  this  purpose, 
the  reflection  of  the  incident  wave  off  a  rigid,  immovable  body  is  referred  to  as  the  scattered  field, 
and  the  incremental  pressure  field  produced  by  motion  of  the  surface  alone  is  called  the  radiated  field. 
The  wave  field  and  scattered  field  together  are  known  as  the  rigid  body  pressure  field. 

The  subaoue  division  used  in  calculating  the  time-delayed  influences  of  surface  pressures  on 
each  other  varies  with  the  distance  between  pressure  zones.  In  the  current  implementation  the  refine¬ 
ments  of  subzone  meshes  are:  1.  Coarse,  2.  Fine,  and  3.  Singular  zone  meshing.  The  coarse  meshing 
is  used  for  two  zones  that  are  remote  from  each  other;  the  fine  mesh  for  two  zones  near  each  other; 
and  the  singular  zone  mesh  for  calculating  the  influence  of  a  pressure  zone  on  itself. 

INCIDENT  WAVE 

Simple  incident  wave  definitions  can  be  given  for  plane  and  spherical  waves.  The  incident  field 
for  a  plane  wave  is  defined  by  specifying  a  time  history  of  the  incident  pressure  on  one  plane  in 
space,  the  coordinates  of  one  point  on  the  plane,  and  the  direction  of  wave  travel  normal  to  the  plane. 
Then  the  pressure  history  at  any  other  point  is  time-delayed  by  the  travel  time  to  the  point  from  the 
defining  plane.  This  relationship  is  expressed  as 

P(XJ)  -  P0(t  -  ( X-I0yT/c ) 

where  p(f,t)  is  pressure  at  point  1  at  time 
p0  is  the  pressure  history  at  Z0  , 

T  is  the  incident  wave  direction, 
c  is  the  wave  velocity. 

For  a  spherical  wave,  the  pressure  history  is  given  for  one  spherical  mathematical  surface 
surrounding  the  source  of  the  wave,  and  the  coordinates  of  the  source  and  the  radius  to  the  defining 
surface.  Because  spherical  waves,  as  they  spread,  attenuate  in  proportion  to  the  inverse  radius,  this 
factor  must  be  included  in  the  incident  pressure  calculations.  The  pressure  expression  is  as  follows: 

p(X,i)  -  ( R/r)pK(t  -  ( r-R)/c ) 

where  r  ■  |  X-X„  | .  the  distance  from  the  pressure  source, 
and  pK  is  the  pressure  history  at  radius  R  from  the  wave  source. 

All  pressures  prior  to  the  arrival  of  the  incident  wave  are  set  to  zero. 

At  present  only  the  plane  wave  opdon  is  installed  in  the  program,  but  the  spherical  wave  option 
is  easily  added. 

NUMERICAL  RESULTS 

The  two  basic  problems  investigated  for  the  three-dimensional  technique  are  a  submerged  rigid 
immovable  sphere  in  an  infinite  fluid  field,  and  an  elastic  spherical  shell  in  the  same  environment. 
Each  is  subject  to  the  same  planar  step  wave.  No  gravity  is  included  in  the  models.  Comparisons  are 
made  of  different  approximations  of  a  step  wave.  The  problem  dimensions  and  properties  are  as  fol¬ 
lows: 

shell  radius-to- thickness  ratio  a/h  *  50 , 
sheil-to- fluid  density  ratio  p,/pw  ■  7.80, 
modulus  ratio  E/pwc 2  -  97, 

Poisson’s  ratio  v  *  0.3. 


The  structural  model  for  both  problems  is  shown  in  Fig.  2.  This  is  a  quarter  section  of  a  sphere 
used  to  model  an  entire  spherical  shell  with  symmetry  conditions  imposed  on  the  two  cutting  planes. 
Both  the  retarded  potential  program  and  the  ADINA  finite  element  program  allow  for  such  symmetry 
conditions.  The  symmetries  were  imposed  to  save  computation  cost  and  can  be  allowed  because  both 
the  test  structure  and  the  test  loadings  are  axisymmetric  with  respect  to  the  loading  direction.  The 
subdivision  into  pressure  zones  and  finite  elements  is  more  refined  in  the  direction  of  wave  travel  in 
order  to  capture  the  more  detailed  behavior  variations  expected  in  that  direction.  The  particular  zone 
subdivisions  used  in  this  study  are  shewn  in  Fig.  3.  A  fine  mesh  is  used  throughout  the  interzone 
inftu^ncg  calculations. 

The  Rigid  Scattering  Problem 

The  retarded  potential  technique  was  first  studied  for  the  rigid  scattering  problem,  independent 
of  the  structural  analysis  code,  to  verify  the  3-dimensional  integration  method  under  shock  loading. 
Pictorially,  this  is  equivalent  to  removing  the  left  side  of  the  problem  illustrated  in  Fig.  1  and  setting 
the  accelerations  identically  to  zero. 

In  choosing  continuous  approximations  to  a  discontinuous  step  wave  loading,  three  considera¬ 
tions  are:  1.  the  approximation  should  accurately  capture  the  shape  and  peak  values  of  the  surface 
pressures  and  structural  response;  2.  the  loading  should  induce  a  minimal  amount  of  numerical  insta¬ 
bility  in  the  computation;  and  3.  the  loading  data  should  be  easy  to  generate.  The  three  loading 
approximations  tested  in  this  study  are  a  linear  ramp,  one  half  of  a  sine  squared  waveform,  and  a 
Fourier  series.  Because  the  sine  squared  approximation  was  found  to  yield  the  same  solution  behavior 
characteristics  as  the  ramp  function,  the  following  discussion  will  be  limited  to  the  ramp  and  Fourier 
series.  The  motivation  for  trying  the  Fourier  series  approximation  is  to  eliminate  higher  frequency 
components  in  the  loading  that  might  initiate  or  aggravate  numerical  instabilities. 

For  each  loading  type,  variation  can  be  made  in  the  steepness  of  the  approximation  to  an  abso¬ 
lute  step.  For  the  ramp  loading,  this  is  a  variation  in  the  number  of  time  steps  to  reach  peak  load, 
and  for  the  Fourier  series,  a  variation  in  the  highest  frequency  present.  There  is  a  tradeoff  between 
solution  accuracy  and  solution  stability;  the  steeper  loadings  give  more  accurate  but  less  stable 
responses.  A  variety  of  loading  steepnesses  were  tried,  and  the  ones  found  to  give  the  best  comprom¬ 
ise  between  accuracy  and  stability  are:  a  ramp  loading  allowed  to  rise  to  maximum  pressure  in  5  time 
steps,  and  a  Fourier  series  approximation  including  frequencies  up  to  600  Hz  in  1  Hz  intervals.  The 
loadings  are  plotted  in  Fig.  4.  Computations  were  carried  out  to  20  non-dimensional  times  (T=*cr /a) 
to  folly  test  the  stability  of  the  responses.  Because  the  response  in  all  cases  reaches  a  steady  state  by 
five  non-dimensional  times,  the  results  of  the  rigid  body  study  are  plotted  for  this  length  of  time  in 
Figs.  6  with  the  location  of  the  plotted  points  shown  in  Fig.  5.  At  one  location  the  response  is  plot¬ 
ted  out  to  a  foil  20  non-dimensional  times  to  show  the  later  stability  behavior  of  the  computation,  typ¬ 
ical  of  all  locations. 

The  standard  of  comparison  used  in  the  rigid  body  study  is  the  eight-term  Legendre  series  solu¬ 
tion  from  Huang  [11].  As  can  be  seen,  both  loadings  give  accurate  pressure  histories  on  the  surface 
at  all  points.  The  ramp  loading  characteristically  gives  sharper  pressure  peaks  than  the  Fourier  series 
loading  and  higher  maximums.  Theoretically,  the  illuminated  half  of  the  sphere  should  experience  a 
peak  of  twice  the  unit  step  pulse.  The  approximate  loadings  do  not  achieve  this  maximum  because 
the  peak  loading  pressure  is  reached  over  a  period  of  time  rather  than  instantaneously,  allowing  the 
scattered  wave  time  to  dissipate.  The  ramp  loading  characteristically  has  a  high  frequency,  low 
amplitude  oscillation  superimposed  on  the  solution.  This  does  not.  however,  lead  to  divergence.  The 
computations  for  neither  case  diverge  over  the  20  non-dimensional  times  for  which  the  test  was  run. 
The  undulating  approximations  to  zero  pressure  that  the  Legendre  series  solution  and  the  Fourier 
loading  response  exhibit  prior  to  the  actual  arrival  of  the  wave  at  a  given  surface  point,  are  artifacts 
of  the  series  approximations.  This  is  particularly  evident  in  the  higher  8  angles  (Figs.  6). 

A  new  stability  consideration  was  discovered  for  the  computations.  It  was  observed  that  the 
time  step  size  and  the  subzone  mesh  size  should  be  appropriately  matched  to  obtain  stable  solutions. 
The  stability  requirement  was  established  in  previous  studies  that  the  distance  from  an  integration 
point  in  one  pressure  zone  to  the  field  point  in  another  zone  should  be  at  least  a  full  time  step  to 
prevent  an  instantaneous  pressure  influence  from  one  zone  to  the  next.  However,  it  is  now  also 


apparent  that  the  time  step  should  not  be  small  relative  to  the  subzone  integration  elements.  This  is  to 
prevent  lumping  of  integration  element  influences  at  widely  separated  time  steps,  leaving  intermediate 
rime  steps  with  no  apparent  influence.  This  possibility  can  arise  because  the  subzone  element  influ¬ 
ence  is  lumped  at  the  center  of  the  element  prior  to  lumping  the  influences  at  the  discretized  retarded 
times.  The  preferred  relationship  is  to  have  a  time  step  transit  distance  somewhat  smaller  than  the 
closest  distance  between  any  pair  of  integration  and  field  points  which  are  in  different  zones,  but  sta¬ 
bility  problems  occur  only  with  gross  disparity  in  sizes.  This  area  warrants  further  study. 

Submerged  Elastic  Structural  Response 

For  the  combined  problem  of  fluid  and  elastic  structural  dynamics,  the  same  spherical  geometry 
and  plane  step  pressure  loading  is  used  as  for  the  rigid  body  problem.  A  thin  elastic  shell  of 
diameter-to-thickness  ratio,  100,  is  modeled  with  9-noded  quadrilateral  shell  finite  elements.  The 
identical  rnesh  is  used  for  the  finite  element  modeling  of  the  structure  as  for  the  discretization  into 
constant  pressure  zones,  in  accordance  with  the  implementation.  The  simulation  of  structural 
response  to  the  plane  step  wave  is  carried  out  for  a  period  of  20  non-dimensional  times  beyond  the 
incident  wave’s  initial  impact  with  the  sphere.  The  same  aon-dimensional  time  step  of  0.05  is  used  as 
for  the  rigid  body  problem.  A  consistent  mass  formulation  is  used  in  the  finite  element  model  and 
has  been  found  essential  to  the  accuracy  of  the  solution.  No  internal  structural  damping  is  included. 

The  resulting  history  of  surface  displacements  and  velocities  are  plotted  in  Figs.  8  and  9.  Fig¬ 
ure  7  illustrates  the  locations  and  orientations  of  these  plotted  values.  The  angle,  9,  of  the  location 
varies  from  zero  at  the  point  of  of  first  wave  impact  to  180  degrees  at  the  point  of  deepest  shadow. 
All  of  the  surface  displacements  plotted  are  the  components  perpendicular  to  the  direction  of  the 
incident  wave  travel,  and  all  the  velocities  are  the  components  in  the  direction  of  the  incident  wave. 
In  these  figures,  the  numerical  solutions  are  compared  to  a  closed-form  series  solution  for  a  step  wave 
on  an  elastic  sphere  obtained  by  separation  of  variables  from  Ref.  [3.3].  The  same  two  loading  cases 
considered  for  the  rigid  body  case  are  compared:  a  five  time  step  ramp  approximation  of  the  step 
wave,  and  a  Fourier  series  approximation  consisting  of  the  frequencies  from  0  to  600  Hz  in  1  Hz 
increments. 

The  results  show  some  general  trends.  The  results  agree  very  well  in  the  shape  of  the  plotted 
response  curves,  and  have  good  agreement  in  amplitudes  at  early  time.  For  both  loading  approxima¬ 
tions,  the  velocities  and  displacements  exhibit  an  increasing  damping  of  response  with  time  compared 
to  the  series  solution.  The  numerical  solutions  begin  to  increasingly  oscillate,  as  well.  The  oscilla¬ 
tions  are  more  noticeable  in  the  velocities  than  the  displacements.  It  is  evident  that  the  velocity  oscil¬ 
lations  are  integrated  out  in  the  displacements. 

A  few  exceptions  to  these  trends  should  be  mentioned.  The  9*0  velocity  of  the  ramp  loading 
shows  an  exaggerated  oscillation  from  the  beginning,  that,  in  fact,  dies  down  noticeably  before  build¬ 
ing  up  again  at  later  time.  The  Fourier  loading  doesn’t  induce  this,  and  its  response  is  rather  damped 
compared  to  the  initial  spike  in  the  series  solution.  At  both  the  point  of  initial  impact  and  in  the  deep 
shadow  region,  the  tendency  to  oscillate  is  much  stronger  than  elsewhere,  which  can  be  attributed  to 
the  higher  frequency  response  of  the  loading-direction  velocities  at  these  points.  As  one  gets  closer  to 
9  *  90® ,  these  velocities  increasingly  represent  only  the  overall  body  translation,  whereas  at  9  =  0 
and  180,  there  is  a  large  component  of  looil  deformation  included. 

In  comparing  the  response  of  the  two  approximate  loadings,  the  results  generally  are  the  same, 
with  certain  exceptions.  First,  there  is  the  already  mentioned  greater  early-time  oscillations  of  the 
ramp  loading  case.  Then,  there  is  a  difference  in  the  pattern  of  later-time  oscillations.  For  the  ramp 
case,  the  oscillations  tend  to  increase  steadily,  whereas  the  Fourier  series  response  exhibits  a  strong 
pulsation  or  beat  in  the  amplitudes.  This  is  particularly  evident  in  the  velocities  at  9  equal  to  0.  162. 
and  180.  Even  at  early  time,  the  ramp  loading  response  has  a  very  fine  oscillation.  While  the  retarded 
potential  integral  solution  for  the  rigid  body  is  stable  for  the  cases  studied  and  the  Newmark  numeri¬ 
cal  integration  used  in  the  finite  element  analysis  is  unconditionally  stable  for  undamped  structural 
vibration,  the  combined  method  exhibits  divergence,  and  numerical  damping  is  required  to  achieve  a 
solution.  Numerical  damping  by  choice  of  Newmark  parameters  ts  successful  in  stabilizing  the  solu¬ 
tion  without  appreciably  harming  accuracy  at  early  time.  The  numerical  instability,  composed  of 


modi  higher  frequencies  than  the  significant  response,  is  damped  out  faster  than  the  latter.  This  is 
because,  for  a  given  frequency  component  of  the  vibration,  the  damping  causes  the  solution  to  decay 
by  a  fixed  fraction  per  cycle,  and  the  higher  frequency  contributions  to  the  vibration  die  out  faster 
than  the  lower  frequency.  However,  at  later  times  the  stabilized  solution  does  show  a  considerable 
reduction  in  amplitude  compared  to  the  closed- form  solution,  likely  caused  by  the  numerical  damping. 

The  ramp  loading  shows  a  greater  tendency  than  the  Fourier  series  for  instability,  given  the 
wM  set  of  Newmark  parameters.  This  tendency  is  quite  nodceabie  in  the  velocities  at  early  time. 
However,  the  Fourier  series,  in  some  cases,  shows  greater  instability  at  later  times.  The  relative 
instabilities  may  depend  on  which  loading  curve  is  changing  more  rapidly  at  a  given  time.  The  ramp 
has  two  sharp  changes  in  direction  initially  and  then  becomes  constant,  while  the  Fourier  series  con¬ 
tinues  to  fluctuate  throughout  its  loading  history.  By  using  greater  numerical  damping  for  the  ramp 
loading  case,  early  time  stability  equal  to  that  of  the  Fourier  series  can  be  obtained  with  only  slight 
reduction  in  the  amplitudes  of  the  response.  For  the  velocity  at  0  =  0,  the  Fourier  series  still  gives  a 
better  match  of  the  analytical  solution  than  does  the  ramp,  which  oscillates  noticeably  about  the 
closed-form  curve.  However,  the  ramp  more  closely  approximates  the  initial  velocity  peak.  For  the 
plotted  data,  the  Newmark  integration  parameters  used  are  5  ■  0.8  and  a  =  0.4225.  There  is  a 
need  for  further  study  on  the  effect  of  mesh  refinement  and  time  step  size  on  solution  accuracy  and 
stability,  as  well  as  why  the  combined  method  is  unstable  at  all.  A  stability  analysis  is  needed  for  the 
Newmark  integration  procedure  that  Includes  the  presence  of  radiation  damping. 

CONCLUSIONS 

This  study  has  demonstrated  a  folly  three-dimensional  implementation  of  a  fluid-structure 
interaction  analysis  procedure  based  on  the  retarded  potential  integral  coupled  to  a  shell  finite  element 
program.  A  continuous  wave  approximation  to  a  discontinuous  step  wave  has  been  found  to  be  suc¬ 
cessful  in  modeling  shock  response  problems.  The  method  uses  techniques  which  may  be  applied  to 
arbitrary  three-dimensional  bodies,  linear  or  non-linear,  whose  fluid-structure  interface  is  a  closed, 
smooth  surface. 

Although  difficulties  have  been  encountered  in  the  tradeoff  between  overdamping  and  oscillatory 
divergence,  the  method  has  been  demonstrated  to  be  quite  effective.  The  early  time  predictions  for 
both  the  rigid  and  elastic  problems  are  accurate  and  have  the  potential  for  improved  accuracy  up  to 
the  limit  of  the  modeling  refinement  to  predict  higher  frequency  responses.  For  the  geometry  and 
loading.'  tested,  the  rigid  body  problem  has  proven  to  be  persistently  stable,  but  for  the  elastic  prob¬ 
lem  there  is  a  numerical  divergence  at  later  time. 

Further  study  in  the  effects  of  finite  element  grid  size  and  configuration  and  subzone  mesh  divi¬ 
sion  on  problem  accuracy  and  stability  is  needed.  This,  as  well  as  an  investigation  of  the  coupling  of 
the  two  computational  processes,  should  improve  the  performance  of  the  method. 
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Fig.  1  —  Computational  Flow 


Fig.  2  —  Mesh  of  pressure  zones/ finite  elements 


Fig.  3a  —  Non-singular  zone  mesh 


Ftg.  3b  —  Non-singular  zone  mesh 
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